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The reduced dynamics for dark and bright soliton chains in the one-dimensional nonlinear 
Schrodinger equation is used to study the behavior of collective compression waves. For appropriate 
conditions, the reduced dynamics derived from perturbation and variational techniques allows to de¬ 
scribe a chain of dark or bright solitons as a chain of effective masses connected by nonlinear springs 
taking the form of a Toda lattice model on the soliton’s positions. In turn, the Toda lattice possesses 
exact solitary travelling compression wave solutions corresponding to travelling compression waves 
in the original soliton chain. We coin the term hypersoliton to describe such solitary waves riding 
on a chain of solitons. We corroborate our analytical results with direct numerical simulations of 
the nonlinear Schrodinger equation. It is observed that in the case of dark soliton chains, the for¬ 
mulated reduction dynamics provides an accurate description for the robust evolution of travelling 
compression waves. In contrast, bright soliton chains do not have such stable propagating solutions 
due to the desynchronization of the mutual phases between consecutive solitons during evolution. 
Finally, as an application to Bose-Einstein condensates trapped in a standard external harmonic 
trapping potential, we study the case of finite dark soliton chain confined at the center of the trap. 

We find that when the central chain is hit by a dark soliton initiated at the edge of the external 
trap, the energy is transferred through the chain as a hypersoliton. When the hypersoliton reaches 
the end of the central chain, a dark soliton is ejected away from the center of the trap and, as it 
returns from its excursion up the trap, hits the central chain producing again a hypersoliton. This 
periodic evolution is the equivalent of the classical Newton’s cradle. 


I. INTRODUCTION 


Bose-Einstein condensation occurs when a dilute gas 
of bosonic atoms is cooled below a critical temperature 
where a considerable fraction of the atoms occupy the 
same quantum state according to Bose-Einstein statis¬ 
tics. Bose-Einstein condensates (BECs)were first theo¬ 
rized by Bose and Einstein in the 1920s [l| but not exper¬ 
imentally realized until 1995 [2|,[3|[for which the 2001 No¬ 
bel Prize in Physics was awarded [J. Typically, rubidium 
or sodium atoms are used and are cooled to nanokelvin 
temperatures using a combination of laser an evaporative 
cooling. The condensate is held in position by a combi¬ 
nation of magnetic and optical traps. For sufficiently 
low temperatures, the mean field dynamics of BECs in 
a quasi-one-dimensional (ID) trap can be accurately de¬ 
scribed by the so-called Gross-Pitaevskii (GP) equation 
that is a variant of the nonlinear Schrodinger (NLS)equa- 
tion incorporating the external trapping potential [3 • By 
appropriately adimensionalizing time, length and energy 
(see Ref. [5[ for details), it is possible to cast the ID GP 
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equation as 

iut =-luxx +g\u\‘^u-\-VMTU, (1) 

where the rescaled condensate wavefunction is given by 
u{x,t), Vmt{x) is the effective ID (magnetic) trapping 
potential confining the BEG and g = ±I indicates 
whether the atoms have an attractive {g = —I) or repul¬ 
sive {g = +1) scattering length. This ID reduction of the 
system is achieved by the so-called cigar-shaped external 
trapping potential for which two transverse directions are 
tightly confining (such that, effectively, only the ground 
state along these direction is possible) while the longitu¬ 
dinal (in our case x) direction is loosely trapped allowing 
for the dynamics of Eq. (HD to evolve along this direction. 

Since the experimental realization of Bose-Einstein 
condensation, the study of this new form of matter has 
been the focus of intensive theoretical and experimental 
efforts BEGs continue to be a testbed for accessing 

quantum mechanics at a macroscopic level allowing for 
direct observation of matter-wave solitons [El, [8[. Under 
strong transverse confinement in two spatial directions, 
a BEG can be rendered effectively quasi-ID (Ef. In this 
case, depending on the sign of the scattering length be¬ 
tween the BEG entities (usually alkali atoms), it is pos¬ 
sible to observe bright liini (for attractive interactions) 
and dark [iilii (for repulsive interactions) solitons. In 
the present work, we are interested in studying the coffee- 
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tive dynamics of chains of these ID solitons and, specif¬ 
ically, the possibility of stable solutions that coherently 
propagate compression waves along the soliton chain. 

Solitons are ubiquitous nonlinear waves that occur in 
a wide range of physical systems such as plasmas, molec¬ 
ular chains, optical fibers, and long water waves M- In 
many physically relevant setups solitons are extremely ro¬ 
bust (with respect to parametric perturbations) and sta¬ 
ble (with respect to configuration perturbations). They 
can interact elastically with other solitons, travel long 
distances, and travel through inhomogeneities with min¬ 
imal deformation and dispersion. This striking stability 
relies on the balance between dispersion and nonlinear¬ 
ity. For instance, in the absence of external trapping 
{Vmt = 0) and for the case of an attractive condensate 
{g = —1), the homogeneous background GP ([B accepts 
exact bright soliton (BS) solutions of the form^, [13 

Mbs = asech(a(a; — ^(t))) e® (2) 

where a is the amplitude (and inverse width) of the soli- 
ton, ^(t) = vt^^o is its position (^o being its initial loca¬ 
tion) and its phase is given by 0bs(^) = (a^ — i;^)t/2 + 

(00 being its initial phase). On the other hand, in the 
case of a repulsive condensate, for a homogeneous and 
stationary background density, the GP equation ffl ac¬ 
cepts exact darA; soliton (DS) solutions of the form 

Wds = \/«0 [B tanh ({x - + iA] (3) 

where no is the density of the constant background that 
supports the DS, ^(t) is its position as defined above for 
a BS and its phase is given by 0ds(A) = —not + 0o (0o 
being its initial phase). The parameters of the DS are 
related by the following expressions: = 1 and 

V = Ay/n^. 

In this work we study the collective dynamics of chains 
of interaction bright and dark solitons. The manuscript 
is organized as follows. The next section is devoted to 
describing the dynamical reduction where the evolution 
of a chain of well separated and nearly identical solitons 
can be reduced, for both BSs and DSs, to a chain of 
effective particles connected with nonlinear springs mod¬ 
elled by the fundamental Toda lattice on the soliton’s 
positions. Section [Till is devoted to constructing appro¬ 
priate initial conditions for the original GP equation to 
support Toda lattice solitons riding on chains of BSs and 
DSs, a.k.a hypersolitons. We present numerical results 
from direct integration of the GP equation which very 
closely match the solutions of the corresponding Toda 
lattice prediction. We describe the robustness of the con¬ 
structed hypersoliton solutions and present some typical 
collision scenarios. Also, in this section, motivated by 
the presence of harmonic trapping in typical BEG exper¬ 
iments, we study the effects of considering a finite soliton 
chain that is confined at the bottom of the external po¬ 
tential. Specifically, we present numerical results for a 
finite chain of DSs supported by an external harmonic 
trap giving rise to dynamics that are akin to the oscilla¬ 
tions of the classical Newton’s cradle. Finally, in Sec. [IVl 


we summarize our results and present possible avenues 
for further research. 


II. DYNAMICAL REDUCTION FOR SOLITON 
CHAINS 

In this section we summarize the dynamical reduction 
for chains of BSs and DSs. Under suitable conditions, 
both systems reduce to a Toda lattice on the position 
of the solitons. Hence, by initializing the soliton train’s 
positions and velocities according to the Toda lattice soli¬ 
ton, a travelling compression pulse can be sustained. 

A. Bright soliton chains 

The BS solution (|2]) describes the coherent evolution 
of a density heap on a zero background in an attractive 
quasi-ID BEG in the absence of an external confining po¬ 
tential. When an external trapping potential is included 
and/or in the presence of other BSs, the BS is perturbed 
inducing a deformation of its shape. However, under 
small perturbations, and noting that BSs are robust, it is 
possible to approximately describe the dynamics of the 
BS by the ansatz m as long as its amplitude, width, 
position, velocity and phase are dynamically adjusted to 
follow accurately the actual solution of the system. For 
instance, in the presence of a magnetic trap of the form 

Vmt{x) = - (4) 

where the strength of the trap D is small (compared to 
the soliton width), a single BS solution will undergo left- 
to-right periodic oscillations of frequency U |2QI423| . On 
the other hand, the presence of another BS, provided 
that both solitons have similar amplitudes and veloci¬ 
ties and that their separation is large compared to their 
widths [ensuring that their shape can still be approxi¬ 
mated by Eq. their interaction dynamics can be re¬ 
duced to a set of coupled ordinary differential equations 
(ODEs) mil. These reduced ODEs depend on all the 
parameters of the solitons. Namely, defining a vector of 
time dependent parameters Pi = (a^, 00 for the 

Tth soliton containing, respectively, its amplitude, posi¬ 
tion, velocity and phase, the dynamics for the BSs can 
be approximately described (under the above mentioned 
conditions) by a set of coupled ODEs on the parameters 
Pi as follows 

CLj = 4:0 — Sj^j 

< = Vj — 2{Sj^j-i-\-Sj^jj^i)^ (5) 

, + 6 Vj (Cj, + Cj, ), 
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where 


c. 

^j,n 

Cj,n 
4^j,n 


^j ^n Cn)? 

1 ~ ~<^i,i+i- 


(6) 


For our particular consideration, we are interested in a 
chain of well-separated, almost identical (a^ ~ a), BSs. 
Under these conditions, the equations of motion for an 
infinite chain of ordered {^i < Ci+i) BSs can be further 
simplified to (see Ref. and references therein) 

= cri -14 4a^ - (7^,^+! 4a^ (7) 


where aij = ±1 is determined by the relative phase of 
consecutive (|i — j\ = 1) BSs. Namely, > 0 corre¬ 
sponds to out-of-phase (OOP) BSs (|0j — 4>i\ = tt) and 
aij < 0 corresponds to in-phase (IP) BSs {(pj — (pi = 0). 
This means that OOP BSs experience mutual repulsion 
while IP BSs experience mutual attraction. It is evident 
that a homogeneous, equidistant, chain of BSs described 
by Eq. © is a fixed point of the system. Furthermore, 
it is straightforward to prove that this fixed point in the 
reduced model is unstable in the case of IP BSs and it 
is neutrally stable for OOP BSs. However, stability of 
the fixed point in the reduced model does not imply sta¬ 
bility of the steady state of the original GP system. In 
fact, stability of the reduced model is a necessary but 
not sufficient condition for stability of the original GP 
system. 

Since we are interested in excitations of the homoge¬ 
neous chain, let us first study their stability. In Fig. [U 
we depict typical time evolutions corresponding to an IP 
(top-left) and OOP (top-right) BS chain. The system 
is initialized using the numerical steady state of the GP 
system with periodic boundary conditions —found using 
a standard fixed point iteration algorithm— with a small 
perturbation. As can be seen in the figure, both the IP 
and OOP cases are unstable. The nature of the instabil¬ 
ity seems, however, different. The IP case, for which we 
know that even in the reduced ODE model is unstable 
due to the mutual attraction between BSs, seems to be 
strongly unstable. The instability is manifested by two 
neighboring BS coalescing as early as t ^ 50. This in¬ 
stability is easy to understand since a small perturbation 
will induce two neighboring BSs to be slightly closer than 
its other neighbors, thus accelerating the process of at¬ 
traction and hence leading to a rapid collision between 
these two BSs. On the other hand, in the OOP case, the 
reduced chain is neutrally stable and, thus, perturbations 
with respect to the positions of the BSs from the equidis¬ 
tant chain should not cause instabilities. As shown in the 
top-right panel of Fig. [U this is the case for intermedi¬ 
ate times {t < 500) where the mutual repulsion between 
BSs is responsible for collision avoidance between neigh¬ 
boring BSs. Nonetheless, as the panel shows, for later 




FIG. 1: Color online. Dynamics of a slightly perturbed 
equidistant chain of 12 BSs of amplitude = 1 in the pe¬ 
riodic domain x E [—54,54], namely a separation of ro = 9. 
The top left and right panels depict, respectively, the spatio- 
temporal evolution of the (square root of the) density \u{x, t)\ 
for the IP and OOP initial chains. The initial conditions are 
set as the numerically exact steady states found by Newton 
iterations and then the BSs positions are perturbed with ran¬ 
domly distributed displacements between —ro/18 and ro/18. 
The dashed lines represents the orbits obtained from the cor¬ 
responding reduced Toda lattice model (|9|) (for the IP case 
this model can only be integrated until the first collision time 
for t < 50). The bottom two panels depict the phase differ¬ 
ence between consecutive density maxima. These correspond 
to the relative phases Ap between consecutive BSs (top and 
bottom subpanels corresponding, respectively, to the IP and 
OOP cases). 


times, t ~ 575, a collision between neighboring BSs does 
indeed occur. The presence of a collision is unequivocal 
evidence that the involved BSs where not OOP when they 
collided. In fact, the loss of the OOP property between 
consecutive BSs is precisely what induces the instability 
of the otherwise OOP initial BS chain. The desynchro¬ 
nization between the phases can be clearly seen in the 
bottom two subpanels of Fig. [TJ The panels depict the 
time evolution of the relative phase between consecutive 
BSs that initially start in the IP (top subpanel) and the 
OOP (bottom subpanel) configurations. It is clear that 
the OOP property between consecutive BSs is approx¬ 
imately held for intermediate times (see t < 500 in the 
bottom subpanel) ensuring the mutual repulsion between 
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FIG. 2: Color online. BdG stability spectrum for the steady 
state for a chain of 10 equidistant chain of IP (top left) and 
OOP (top right) BSs in the periodic domain x G [—35,35]. 
The bottom set of panels depict, from top to bottom, the 
steady state and the first nine most unstable eigenfunctions 
where the real part is depicted by the (blue) solid line and 
the imaginary part by the (red) dashed line. 


consecutive BSs and, thus, stability for the chain. How¬ 
ever, it is also clear that the OOP property gets pro¬ 
gressively worse until a pair of consecutive BSs have a 
zero relative phase between them, i.e. they are IP around 
t = 575, inexorably leading to their collision. 

To further understand the nature of the instability, we 
have computed the Bogolyubov-de Gennes (BdG) stabil¬ 
ity spectrum for the IP and OOP steady states of the 
GP system. The BdG spectrum is computed by consid¬ 
ering perturbations from the steady state uo{x) of the 
following form: 

u{x,t) = ^i4o(a:^) -\-e{a{x)e^^ + 6(x)e“^ ^)^ , (8) 

where fi is the so-called chemical potential (or [the neg¬ 
ative of] the temporal frequency) of the solution, A is 
the eigenvalue with associated eigenvector {a, 5*}, and 
where (•)* stands for complex conjugation. After apply¬ 
ing the perturbation ansatz (|8j) into the GP equation 
(HD and linearizing the ensuing equation, one obtains 
an eigenvalue problem with a corresponding eigenfunc¬ 
tion w{x) = a(x) + 6*(x) at t = 0. After computing 


the spectrum, any eigenvalue A = + i\i with a pos¬ 

itive real part (A^) indicates an unstable eigenfunction. 
The spectrum associated with the IP and OOP steady 
states is depicted, respectively, in the left and right top 
panels of Fig. O As is can be noticed, the OOP spec¬ 
trum has a handful of complex eigenvalues with a small 
{Xr < 0.02) real part indicating a weak instability. In 
contrast, the IP case reveals a larger (purely real) in¬ 
stability with max(A^) ^0.123, indicating a stronger in¬ 
stability. Closer inspection of the unstable eigenfunctions 
(see bottom set of panels in Fig.[2j) reveals that the insta¬ 
bilities manifest themselves as local translational modes 
for consecutive BSs in the opposite direction and, thus, 
bringing them closer to each other. We have checked 
that the stability results above (cf. Figs. [T] and [2]) are 
very similar for other values of the parameters such as 
amplitude, number, and separation of the BSs, as well 
as different domain lengths. Evidently, as more BSs are 
included in the system, a higher degeneracy of the eigen¬ 
values arises since all solutions and eigenfunctions posses 
translational symmetry. Finally, it is relevant to mention 
that the BS chain might be rendered stable by a suitable 
choice of periodic lattice potential providing stabilizing 
pinning for each BS located at the respective minimum of 
the lattice potential Bin. However, we do not explore 
this avenue further in this manuscript. 

Since the IP BS chain is highly unstable, we will focus 
our attention in the case of the weakly unstable OOP 
chain. Therefore, let us consider the case aij = +1 for 
which the reduced equation of motion yields 

Ci = (9) 

which has precisely the form of the celebrated Toda lat¬ 
tice [32I that is further described below. It is important 
to mention at this stage that the restriction of locked 
phases between BSs (IP or OOP) is not necessary to ob¬ 
tain a Toda lattice-type model. In fact, by allowing the 
phase of each BS to dynamically evolve, the equations of 
motion (|5]) reduce to a complex Toda lattice of the form 

m- 

(ij = 2a , ( 10 ) 

where the corresponding complex variable qj for this 
complex Toda lattice is defined through the original BS’s 
parameters by 


1j = « Ci - j In (a^) + « {j TT - + Sj + S) , (11) 

where a and v are the ensemble average height and ve¬ 
locity of the BSs, while Q is the position of the j-th BS 
and the 6j^s are the BS phases and 6 is their average. 


B. Dark soliton chains 

As was done in the previous section for the BS chain, 
the case of DS chains can also be reduced to a set of 



































5 


ODEs on the DS parameters. This reduction is obtained 
through DS perturbation theory as described in detail 
in Refs. [H, [iqI, In particular, considering a DS 

ansatz of the form m for all DSs in a chain supported 
on a constant background with density no, the equations 
of motion are approximately reduced to 

_ g^3/2 ^ ^^ 2 ) 

which, as for the BS chain, is a form of the Toda lat¬ 
tice [33 on the DS positions. 

The main difference between the reduced dynamics of 
bright and dark solitons that is worth pointing at this 
stage is that BSs can have mutual interactions that are 
repulsive or attractive depending on their relative phases 
as described above. On the other hand, DSs are al¬ 
ways repulsive. Therefore, the stationary homogeneous, 
equidistant, DS chain is always stable. This observation 
will be crucial when comparing the dynamics of BS and 
DS chains from the GP model m and their respective 
dynamical reductions © and m- 


C. Validation of the Toda lattice reduction 

We now validate the dynamics obtained from the 
dynamical reduction for both bright and dark soliton 
chains. In order to numerically approximate an infi¬ 
nite lattice, we take a periodic domain in the inter¬ 
val X G [—L, L] and place the solitons equidistant from 
each other accounting for the periodic boundary condi¬ 
tions. Thus, considering N solitons gives a steady state 
equidistant configuration such that — ^i\ = vq for 

i = 1,...,V — 1, where the distance | • | is measured in 
the periodic domain such that I^at — -fil =2L — r^. 

Let us first consider a BS chain. As we described ear¬ 
lier, the perturbed equidistant chain evolves as depicted 
in Fig. [ 1 ] where the top left and right panels correspond, 
respectively, to the IP and OOP cases. We have superim¬ 
posed the corresponding dynamics of the reduced ODE 
model of the BSs ([9]) using dashed lines. As can be ob¬ 
served from the figure, in the IP case (top left panel) 
before the collision of the BSs (t < 50), the ODE model 
very closely follows the BS dynamics. After collision, the 
ODE model breaks down as the BS centers coalesce and, 
thus, we only show the reduced ODE orbit up to the 
first collision time. On the other hand, for intermediate 
times, the OOP chain (top right panel) does not suffer 
from the collision of BSs as it assumes that all BSs are al¬ 
ways OOP. The resulting reduced ODE dynamics closely 
follows the original GP dynamics for short times, but 
later deteriorates since, as explained earlier, the BSs lose 
the OOP synchronization. Nonetheless, for intermedi¬ 
ate times, while the BSs are kept separated, the reduced 
ODE model does reproduce the original GP dynamics. 

In contrast to the BS chain, the DS chain does not 
suffer from phase-induced instabilities since DSs always 
repel each other. As a result, the reduced ODE model 
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FIG. 3: Color online. Interaction dynamics for 12 DSs ini¬ 
tialized at random positions with random velocities in a peri¬ 
odic domain. Depicted is the spatio-temporal evolution of the 
(square root of the) density. The dashed lines represent the 
orbits obtained from the corresponding reduced Toda lattice 
model ([T2|. 


for the DS chain m provides a very robust model for 
the GP dynamics under extended time evolutions for 
any initial condition provided the DSs are initially well- 
separated. An example of the dynamics from the reduced 
Toda lattice ODE ([T2j) and the original GP model is de¬ 
picted in Fig. [31 where a collection of DSs placed at ran¬ 
dom locations with random initial velocities evolves in 
time. It is clear from the figure that the Toda lattice 
model gives an accurate prediction of the DS positions 
for the original GP system for long times. 


III. TODA LATTICE SOLITONS 

A. Preliminaries 

Before constructing Toda lattice solitons on the chains 
of bright and dark solitons, let us review the form of these 
solutions for completeness. The Toda lattice is one of the 
most popular models in physics since, by construction, it 
was designed as to prescribe a chain of nonlinear oscilla¬ 
tors with completely integrable evolution [^. As such, 
the Toda lattice possesses some exact solutions that are 
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the foundation for building more complex solutions. In 
particular, the Toda lattice possesses periodic and local¬ 
ized solutions [32|. Here we focus on the latter type of 
solutions referred to as Toda solitons. The Toda lattice’s 
equations of motion 


Vn = hTL(^n+l “ Vn) “ hTL(^n “ Vn-l), 
— ^^-biVn-Vn-l) _ ^^-biVn + l-Vn) ^ 


(13) 


of soliton amplitude a for the BS chain and in terms of the 
background density no for the DS chain. It is worth men¬ 
tioning that the uniform pre-compression experienced by 
the periodic chain effectively corresponds to a rescaling 
on the strength of the Toda lattice potential A. This is 
evident when rescaling the soliton positions by a factor 
7 , then the exponential interaction terms be¬ 
come H = ie-Kyr.-yr.-i) 

where A = Ae~^^. 


originate from the interaction of nearest neighbors in a 
one-dimensional chain of coupled, unit mass, particles at 
positions interacting through the potential 

Vrhi^y) = ^ ^y- ( 14 ) 

Here, Ay is the separation between particles and A and 
b are positive parameters prescribing, respectively, the 
strength and decay of the inter-particle interactions. By 
following the evolution of the particles through their mu¬ 
tual separation 


Ayn — ^n+l Uni (1^) 

and defining Sn = d{Ay^)/dt and Pn = dyn/dt so that 
Sn-i — Sn = Pn^ the equations of motion can be rewritten 
in term of 5'^ = f Sndt as 


In 


1 + ^ 
a 


— (5'n+l ~ 25'^ + Sn-1 

m 


). 


(16) 


Then, it is straightforward to find solitary kink solutions 
for this system in the form: 


= tanh {nn ± /3t) + const. 


(17) 


where the kink velocity is c = jd/n and its amplitude p 
is given by 


yd = a/] 46 sinh/D, (18) 

where the width of the kink is a free parameter. It 
should be noticed that this solution is stable and it cor¬ 
responds to a compression wave that travels through the 
lattice js^ . 


B. Toda lattice solitons: hypersolitons 

We now seek to use the soliton solution for the Toda 
lattice (see previous section) to construct a Toda lattice 
soliton on the reduced lattice equations for the bright and 
dark soliton chains. Let us consider an equilibrium con¬ 
figuration consisting of a chain of N equidistant solitons 
with separation tq = 2L/N in the periodic domain x G 
[—I/,L]. Both OOP[i9| BS and DS chains are reduced, 
respectively, to the Toda lattice chains ([9j) and ([12]) where 
the Toda lattice potential parameters are given in terms 


Let us start by initializing the chain of BSs such that 
the initial positions and initial velocities satisfy the cor¬ 
responding Toda soliton (ED- An example of this case is 
depicted in Fig. |4]for = 10 BSs in the periodic chain 
X G [—43,43]. The top panels depict the initial condition 
for the displacements from equilibrium between solitons 
= Cn — Cn-i “ (^0 (left subpancl) and their respective 
initial velocities r^. As it is evident from these pan¬ 
els, the kink dm corresponds to a localized compression 
wave for the soliton’s positions. It is worth mentioning at 
this stage that contrary to the homogeneous, equidistant, 
chain where 2L = N vq, for the chain initialized with the 
Toda lattice solitons the length of the domain has to be 
adjusted since we are introducing a compression wave to 
the initial condition. Thus, we compute the length of the 
domain by adding all the separations between consecu¬ 
tive solitons. The bottom row of panels in Fig. 0] depicts 
the evolution of the density after seeding the original GP 
equation with the Toda lattice soliton initial conditions 
on a pre-compressed BS chain. As it can be observed 
from the left subpanel, the initialized compression wave 
travels at the prescribed speed —for this panel we set the 
final time to precisely 2c/(21/), namely, the time needed 
to perform exactly two complete cycles through the lat¬ 
tice. In the figure, the dashed lines correspond to the 
Toda lattice soliton from the reduced ODE model (|9|) 
which, for short times, accurately approximates the full 
GP dynamics. However, as we have described before, 
the BSs will eventually lose their OOP synchronization 
leading to the coalescence of two consecutive solitons. 
The first coalescence occurs approximately at t = 125 
for the right-most two solitons of the lattice. Nonethe¬ 
less, after this coalescence, the Toda lattice soliton seems 
to reform again which, in turn, suffers from the coales¬ 
cence of more BSs as time progresses. It is interesting 
to note that although the reduced ODE should fail after 
a BS pair loses its OOP synchronization, and as it was 
noted in Ref. for OOP initial conditions, the reduced 
ODE model still captures the dynamics of the interact¬ 
ing chain for some time. Nonetheless, as can be seen in 
the bottom right panel of Eig. 01 after longer integration 
times (10 full cycles around the lattice), the Toda soliton 
does not preserve its shape and other excitations start 
populating the dynamics, including Toda solitons that 
apparently move in the opposite direction of the original 
Toda soliton. By the same token, it is also evident that 
the ODE description (see dashed line) fails to capture the 
long term dynamics of the original GP system. This is 
evidence that the unstable character of the BS chain due 
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FIG. 4: Color online. Toda lattice soliton riding on a chain 
of = 10 BSs in the periodic domain x E [—43,43]. The 
top panels depict the initial conditions in terms of the rela¬ 
tive displacements from equilibrium rn{t = 0) (left subpanel) 
and their corresponding velocities rn{t = 0) (right subpanel). 
The bottom left and right panels depict the evolution of the 
(square root of the) density of the original GP model initial¬ 
ized with the Toda lattice soliton for, respectively, a total time 
equivalent to two and ten cycles of the Toda lattice soliton 
around the periodic domain. The dashed lines correspond to 
the reduced ODE model of the Toda lattice 


to the desynchronization of the phases precludes satis¬ 
factory modeling of the full GP system with the reduced 
ODE (|9|) where all BSs are assumed to be OOP. 

Let us now turn our attention to the DS case. Here, 
the problem of the desynchronization of the phases is nat¬ 
urally avoided and thus it should be possible to observe 
Toda solitons travelling stably through the DS chain. In 
fact, this is precisely what we observe in our numerics 
for a wide range of parameters for the DSs themselves 
and of the Toda soliton. Figure [5] depicts three of such 
examples corresponding to three difference Toda soliton 
initial velocities. The different velocities are computed 
from the Toda lattice soliton solution (p!7|) corresponding 
to the reduced Toda model for the DSs m for no = 1 
and unperturbed initial separations vq = 6, 5, and 4. 
The top panels in the figure depict the initial compres¬ 
sion wave (vn) and its initial velocity (r^). The mid¬ 
dle row of panels depicts the dynamics on the (square 
root of the) density (|n|) together with the reduced ODE 
Toda model ([T^ superimposed to it. The bottom row of 
panels depicts the time derivative of the (square root of 
the) density {\u\t) illustrating that outside of the region 
of localization of the Toda soliton, there are no pertur¬ 
bations indicating a stable, radiationless, propagation of 
the Toda soliton. From now on, we dub such a solu¬ 
tion a hypersoliton as it is a (Toda lattice) soliton riding 



FIG. 5: Golor online. Toda lattice soliton riding in a chain 
of = 12 DSs on a periodic domain. The top two rows of 
panels depict the initial conditions in terms of the relative 
displacement from equilibrium of the mutual distances rn(0) 
(top subpanels) and their corresponding velocities (bottom 
subpanels). The middle and bottom row of panels depict, 
respectively, the evolution of the (square root of the) density 
\u{x,t)\ and its time derivative for the original GP model 
initialized with the Toda lattice soliton. The left, middle and 
right columns correspond to increasingly large Toda lattice 
velocities for a background density no = 1 and, from left 
to right, ro = 6, ro = 5, and ro = 4, which correspond, 
respectively, to c = 0.0180, c = 0.0489, and c = 0.1329. The 
dashed lines correspond to the reduced ODE model of the 
Toda lattice m- 


on a chain of (dark) solitons. As shown in Fig. [5l the 
reduced Toda lattice and, in particular, its Toda lattice 
solution, accurately describes the behavior of the original 
GP model. 

In order to further study the stability and robustness 
of the crafted hyper solitons, we proceed to initialize the 
lattice with an initial condition corresponding to a per¬ 
turbed hypersoliton. Specifically, Fig. [6] shows the dy¬ 
namical evolution after using 30%, 50%, and 80% ran¬ 
dom initial perturbations to, both, initial relative po¬ 
sitions and their respective velocities. It is remarkable 
that the addition of such high levels of perturbations 
to the initial conditions do not destroy the hypersoli- 
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FIG. 6: Color online. Stability of the hypersoliton. The 
layout of the different panels is the same as in Fig. [5] The 
initial condition corresponds to the case of the middle column 
of panels of Fig.[5](ro = 5) with random perturbations added 
to the initial condition. The initial condition is perturbed by 
adding a random perturbation of size, from left to right, 30%, 
50%, and 80% on the initial separations and velocities of the 
DSs in the chain. 


ton. Instead, these perturbation seem to just provide 
some background noise over which the hypersoliton rides 
with minimal interaction. This background noise is more 
clearly visible in the bottom row of panels depicting the 
time derivative of the (square root of the) density. This 
remarkable robustness of the hypersoliton even after the 


addition of such large perturbations to the initial con¬ 
ditions —that in turn develop into perturbations along 
the whole domain— is due to two independent facts: (a) 
on the one hand, DSs for the GP equation are very ro¬ 
bust, a property owing from their topological charge, and 
(b) the structural stability of Toda solitons of the Toda 
lattice. These two stability properties, at two distinct 
levels of the model, combine to give the hypersoliton on 
the original GP model its robustness. 

C. Toda lattice solitons: hypersoliton collisions 

Now that we have established the existence and stabil¬ 
ity of the hypersoliton solutions in DS chains, it is pos¬ 
sible to study several dynamical aspects arising at this 
higher structural level. For example, one can initialize 
the DS chain with two (or more) hyper solitons and al¬ 
low them to interact. It is expected that the dynamics 
governing the interactions between hypersolitons will be 
prescribed by the corresponding dynamics on the Toda 
lattice. In particular, Toda lattice solitons collide elas¬ 
tically without energy loss during the collision process. 
This is precisely what we observe for a wide range of 
cases when seeding two hypersolitons at different loca¬ 
tions with different initial speeds (if they have the same 
velocity they will chase each other indefinitely). Figure [7| 
shows typical examples of hypersoliton collisions. Specif¬ 
ically, the left, middle and right cases correspond, respec¬ 
tively to: (a) a head-on collision for opposite velocities 
with the same magnitude, (b) a head-on collision for ve¬ 
locities with different magnitudes, and (c) two chasing 
hypersolitons where a faster one chases a slower one until 
they collide. As the bottom panels show, all these colli¬ 
sions, as expected, are elastic (i.e., no energy is lost from 
the travelling hypersolitons before or after the collisions). 
It is interesting to note, as per the “particle” nature of 
the hypersoliton structures, owed to their nonlinear char¬ 
acter, there is a shift in their paths when comparing their 
straight line trajectories [in the (x, t) plane] before and 
after collision. This effect is clearly visible in the last 
case of Fig.[7]where the fast soliton is advanced after col¬ 
lision while the slow one is delayed after collision. This 
effect will be important as we construct a quantum ana¬ 
log of the classical Newton’s cradle using hypersolitons 
in a chain of DS trapped inside a confining potential in 
the next section. 


D. Quantum Newton’s cradle 

We now proceed to study the effects of adding a con¬ 
fining external potential [i.e., Vmt 7 ^ 0 in Eq. O], rele¬ 
vant to the modelling of magnetically trapped BECs, on 
the dynamics of hypersolitons supported by a finite DS 
chain. In the presence of the external parabolic potential 
m, a single DS exhibits approximately harmonic oscilla¬ 
tions with a frequency uj = Q/a /2 (see Refs. 
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FIG. 7: Color online. Hypersoliton collisions. The left, mid¬ 
dle and right columns depict, respectively, the collision for (a) 
a head-on collision with equal but opposite velocities, (b) a 
head-on collision with different speeds, and (c) a chasing col¬ 
lision where a faster hypersoliton (seeded on the left) chases 
a slower hypersoliton (seeded on the right) until they collide 
in a nonlinear fashion. All panels have the same meaning and 
layout as previous figures. 


[ 3 ^ and references therein). This result is valid in the 
so-called Thomas-Fermi limit corresponding to the high 
density limit. Therefore, combining —through perturba¬ 
tion theory— the mutual interactions between DSs and 


the force exerted by the external trap yields [110 

(19) 

corresponding to the Toda lattice described by Eq. m 
with an added on-site force generated by the external 
potential on each of the DSs of the chain. The above 
model has not only been validated numerically, but it has 
been used to predict the normal modes of vibration for 
a small number of DSs in actual experiments [^, |40|. In 
fact, the Toda lattice with the on-site potential m pos¬ 
sesses a steady state solution emerging from the balance 
of the mutual repulsions between the DSs and the attrac¬ 
tion of the external trap towards the trap’s center. This 
compressed steady state train of N DSs has N distinct 
normal modes of vibration corresponding to the normal 
modes of vibration of N coupled masses through (linear) 
springs with springs at each end attached to rigid walls. 
For example, for N = 2 there exist 2 normal modes of 
vibration corresponding to the in-phase and out-of-phase 
modes of vibration of the DSs |^ . 

Instead of studying further the normal modes men¬ 
tioned above, we opt here to emulate the dynamics of a 
classical Newton’s cradle using a chain of solitons. The 
idea is to start with an initial stationary chain of N DSs 
at the bottom of the parabolic trap and then drop a sin¬ 
gle, outer, DS from a position higher up in the trapping 
potential. This outer DS will experience the force of the 
trapping potential and ride down the external trap to 
collide with the stationary DS chain. The collision ex¬ 
cites a moving hypersoliton within the inner DS chain. 
When the hypersoliton reaches the opposite extreme of 
the inner DS chain, a single DS is expelled outward. This 
new outer DS will ride up and down the trap until it 
hits the inner DS chain thus repeating the process of a 
soliton analogue of the classical Newton’s cradle. It is 
relevant to note at this stage that a similar idea was ex¬ 
perimentally achieved in a ID Bose gas of ^^Rb atoms 
by initially splitting a wavepacket into two wavepackets 
with opposite velocities. These packets then evolve by 
going up and down the trap and periodically colliding at 
the center 0. Also, in Ref. [42|, the authors propose 
another quantum analogue of a Newton’s cradle using 
a BEG by partitioning the wavefunction using a peri¬ 
odic optical lattice potential. In contrast, our method to 
create a Newton’s cradle is based on the effective non¬ 
linearity of the GP model and the repulsive interaction 
dynamics between dark solitons. 

Eigure [8] depicts three different attempts at recreat¬ 
ing a Newton’s cradle-type evolution with our setup. In 
these examples, we use a stationary DS chain of = 12 
DSs placed at the center of a magnetic trap of strength 
Q = 0.01, namely uj = 0.01/\/2 [see Eq. The sta¬ 

tionary inner DS chain is obtained by a fixed-point iter¬ 
ation (Newton’s) method initialized with a chain of DSs 
with zero velocities positioned at the steady state loca¬ 
tions. Once the stationary inner lattice is found, an ex¬ 
tra, outer, DS is seeded away from it. The distance of the 
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FIG. 8: Color online. Hypersoliton Newton’s cradle. A DS 
soliton (leftmost one) is released at various distances away 
from a stationary lattice of 12 DSs placed at the center of 
the parabolic trapping potential (|3|) with Q = 0.01. The DS 
is released at a distance equivalent to (a) 3ro (left column 
of panels), (b) 12ro (middle column of panels), and (c) 20ro 
(right column of panels) where ro is the separation between 
the central DSs. All panels have the same meaning and layout 
as previous figures. 


outer DS to the inner DS chain is varied and the resulting 
evolution is analyzed. The three cases depicted in Fig. [8] 
correspond to, from left to right, an initial distance of the 
outer DS of (a) 3ro, (b) 12ro, and (c) 20ro where ro is 
the distance between the two innermost DSs of the cen¬ 
tral chain. For case (a), corresponding to a short drop¬ 
ping distance of the DS, the Newton’s cradle dynamics 
is observed for a couple of periods but apparently the 
outer DS is “absorbed” by the inner chain resulting in a 
larger inner chain (i.e., N -\-l DSs) that simply oscillates 
in the in-phase normal mode. The mechanism whereby 



FIG. 9: Golor online. Schematic of the mechanism for the loss 
of the hypersoliton Newton’s cradle through the synchroniza¬ 
tion of the outer DS (red) with the central DS chain (black), 
(a) Release of a DS from the left with a stationary DS chain in 
the central region, (b) After the hypersoliton travels through 
the chain, the rightmost DS is ejected to the right. The re¬ 
maining central DS chain is now asymmetric with respect to 
the center (see vertical dashed line) and starts moving to the 
right, (c) Both outer DS and central DS chain perform half 
an oscillation after their respective excursion up the trap, (d) 
After the DS collides and transfers its energy (through a prop¬ 
agating hypersoliton through the central DS chain) the left¬ 
most DS is ejected. The central DS chain still has the extra 
energy of moving towards the left, (e) Both the DS and the 
central chain now oscillate in the same direction (although 
slightly out-of-phase). This process continues until the en¬ 
ergy of the outer DS is completely transferred to the central 
chain and all that remains is a central chain undergoing left- 
to-right oscillations in the in-phase normal mode. 


the outer DS is absorbed by the inner lattice hinges on 
the fact that, although DS collision are elastic, there is 
a shift in the path of the DSs with respect to the before 
and after collision trajectories as was shown before (see 
discussion on the last collision depicted in Fig. [7]). The 
details of the outer DS “absorption”, or rather the energy 
exchange between the outer DS dynamics and the inner 
DS chain, is explained in Fig. [9l This energy transfer 
is more clearly visible in the left-bottom panel of Fig. [8] 
depicting the time derivative of the (square root of the) 
density. In this panel, it is clear that during the first hy¬ 
persoliton excursion through the inner DS chain, there 
is practically no scattering of energy in the inner chain. 
However, as explained in Fig. [9l just after the hypersoli¬ 
ton ejects the first DS, the inner chain has an extra DS 
on the left and is missing one DS on the right and thus, 
it is no longer close to equilibrium and starts to oscil¬ 
late. This is clearly visible in the bottom-left panel of 
Fig. [8] for 300 < t < 500 where the inner chain moves 
synchronously. This energy transfer mechanism contin- 




















11 


4 



FIG. 10: Color online. Long term evolution of the hypersoli- 
ton Newton’s cradle corresponding to the case depicted on 
the middle column of Fig. [8l Notice the beating of energy 
exchange between the outer DS and the central chain. All 
panels have the same meaning and layout as previous figures. 

ues until all the energy of the outer DS is completely 
depleted and the outer DS is absorbed by the inner lat¬ 
tice that oscillates with its in-phase normal mode (results 
not shown here). 

In order to avoid, or minimize, the energy transfer be¬ 
tween the outer DS and the inner DS chain, it is nec¬ 
essary to decouple the dynamics of the outer DS and 
the inner DS chain. This is achieved by increasing the 
dropping distance of the outer DS. As can be seen in the 
case depicted in the middle column of panels in Fig. [8l 
the excursion of the outer DS has very little effect on 
the inner chain motion. This can be explained because 
in this case, the outer DS travels faster through the in¬ 
ner chain and thus, the latter has less time to develop 
its in-phase normal mode. Furthermore, as the period 
of the outer DS is different from the period of the in¬ 
ner chain in-phase mode, subsequent excursions of the 
outer DS do not synchronize with the in-phase mode. 
The outer DS has a different period in the presence of 
the inner DS chain because its trajectory in (x, t) is the 


concatenation of sinusoids (when the outer DS is trav¬ 
eling up and down the external trap on its own) and 
straight paths (when the hypersoliton traverses the in¬ 
ner chain). [5Q| Therefore, let us effectively consider the 
two involved dynamics, namely the outer soliton oscil¬ 
lations and the inner chain oscillations, as two coupled 
oscillators. These oscillators can synchronize provided 
that their periods are close to each other and that the 
coupling is sufficiently strong • This is precisely what 
happens for a relatively small dropping distance as de¬ 
picted in the first case of Fig. [8l However, as we increase 
the dropping distance, the two coupled oscillators have 
increasingly different frequencies and, at the same time, 
the coupling is reduced as the interaction time between 
the two, given by the time it takes for the hypersoliton 
to traverse the inner chain, is reduced because of a faster 
hypersoliton speed. Thus, in principle, there should be a 
threshold drop-off distance for which the Newton’s cra¬ 
dle should be self-maintained. This is what seems to be 
occurring in the case depicted in the middle column of 
Fig. [8] where apparently a very small amount of energy is 
transferred to the inner chain. To ensure that this small 
transfer does not destroy the Newton’s cradle, we depict 
in Fig. [T0]the same case as in the middle column of Fig. [51 
but for a longer time. As can be seen from the figure, 
there is indeed a transfer of energy from the outer DS to 
the inner chain for t < 9, 000. However, the roles are in¬ 
verted after this time and then the inner chain transfers 
back the energy to the outer DS! This produces a beating 
phenomenon common for coupled oscillators with differ¬ 
ent periods. This periodic energy transfer between the 
outer DS and the inner chain seems to persists for very 
long times (results not shown here) providing a mecha¬ 
nism for a stable, long-lived, Newton’s cradle dynamics. 

Finally, the last column of Fig. [5] depicts an example 
with an even larger drop-off distance. In this case, the 
outer DS also interacts with the edge of the BEG cloud 
and produces sounds waves that remove energy from the 
former and thus, finally settling to a Newton’s cradle 
with slightly lower oscillation amplitude with some back¬ 
ground radiation (sound waves) prevailing in the conden¬ 
sate (results not shown here). 


IV. CONCLUSIONS AND OUTLOOK 

We construct coherent structures consisting of com¬ 
pression waves riding on chains of bright (BS) and dark 
(DS) solitons of the GP model. Namely, a soliton riding 
on a chain of solitons and thus dubbed a hypersoliton. 
The dynamics for chains of BSs and DSs is reduced to a 
Toda lattice on the solitons’ positions, i.e., the solitons 
are modelled as a chain of nonlinearly coupled masses. 
Then, the corresponding Toda lattice solitons (compres¬ 
sion waves on the lattice) can be initialized on the orig¬ 
inal GP model using the exact Toda lattice soliton so¬ 
lution. We show how BS chains are inherently unstable 
due to phase desynchronization between consecutive BSs 
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and thus are poor candidates for supporting hypersoli- 
tons. In contrast, DS chains are stable and DSs, being 
topologically charged, never lose their phases and thus 
are always mutually repelled from each other. We suc¬ 
cessfully craft hypersoliton solutions riding on DS chains 
of the original GP model for a wide range of parame¬ 
ter values. These hypersolitons are robust and stably 
travel at a constant speed without deformation nor radi¬ 
ation. Additionally, we construct multiple hypersolitons 
and observe their elastic collisions in different head-on 
and chasing collisions scenarios. Finally, inspired by the 
classical Newton’s cradle, we study the dynamics of finite 
DS chains trapped inside the customary parabolic exter¬ 
nal potential relevant in experimental Bose-Einstein con¬ 
densates. This is achieved by letting a free, outer, soliton 
hit a stationary inner DS chain creating a hypersoliton 
wave travelling through the latter. As the hypersoliton 
reaches the end of the inner chain, a single DS is expelled 
and allowed to rise and fall down the external trap hitting 
the inner chain repeating the process in a manner akin to 
the classical Newton’s cradle. We study the effects of the 
drop-off distance on the formation of the Newton’s cra¬ 
dle dynamics and argue, in terms of the theory of coupled 
oscillators, that a minimal drop-off distance is required 
for the creation of self-sustained Newton’s cradle oscilla¬ 
tions. 

The present work could be extended in a few inter¬ 
esting directions. For example the effects of finite tem¬ 
peratures in a condensate give rise to dissipation due 
to coupling with the thermal (non condensed) fraction. 


This dissipation can be modelled at the level of the 
GP equation by the so-called phenomenological dissipa¬ 
tion [^, and it is responsible for anti-damping terms 
in the reduced equations of motion of the DSs. It would 
be interesting to analyze the effects of such a dissipative 
term on the dynamics of hypersolitons. On the other 
hand, condensates can be supported by two or more cou¬ 
pled components with linear and/or nonlinear coupling 
terms between them [5[. These coupled models give rise 
to coupled complexes with dark or bright solitons coupled 
to dark or bright solitons in the other component (s), thus 
giving rise to the so called dark-dark and dark-bright 
solitons The dynamically reduced models for 

these coupled systems take the form of coupled Toda lat¬ 
tices [isl . It would be interesting to explore the possibil¬ 
ity to construct hypersolitons and study their stability in 
systems with several components. 
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